Survival of the Fittest: Cox Regression
Capstone 6940
Intro to Cox Regression
Cox regression is a popular regression modeling method used to predict survival rates given certain covariates.
“Survival” can refer to the development of a symptom, time to relapse after remission, or as a time to death [1].
Cox regression model is based on the hazard function h_x(t) with covariates at time t given by:
h_x(t)=h_0(t)\exp(\beta_1x_1 +\beta_2x_2 + \dots + \beta_p x_n) h_x(t) = hazard \ function\ h_0(t) = Baseline \ hazard \ function\
\beta_1x_1 +\beta_2x_2 + \dots + \beta_p x_n = linear\ combination\ of \ covariates\ and\ their\ coefficients
Analysis: Corporate Survival
By Kayla Boyd
Data Source: Employee Turnover Data from Kaggle
| Variable | Description |
|---|---|
| stag | Experience (time) |
| event | Employee turnover |
| gender | Employee’s gender |
| age | Employee’s age (year) |
| industry | Employee’s industry |
| profession | Employee’s profession |
| traffic | From what pipelene employee came to the company. |
| coach | Presence of a coach (training) on probation |
| head_gender | Head (supervisor) gender |
| greywage | The salary does not seem to the tax authorities. |
| way | Employee’s way of transportation |
| extraversion | Extroversion score |
| independ | Independence score |
| selfcontrol | Self control score |
| anxiety | Anxiety score |
| novator | Innovator score |
{r, warning=FALSE, echo=F, message=FALSE}
library(readr)
turnover <- read_csv("turnover.csv")
View(turnover)
library(tidyverse)
library(dplyr)
library(corrplot)
library(RColorBrewer)
library(ggfortify)
library(riskRegression)
library(survival)
library(stringr)
library(zoo)
library(ranger)
library(ggplot2)
library(readxl)
library(MASS)
library(ADGofTest)
library(survminer)
library(car)
# Check for null and duplicate values
cat("Number of missing values :", sum(is.na(turnover)))
cat("Number of duplicates: ", turnover %>%
duplicated() %>%
sum())
# Remove duplicate values
turnover <- unique(turnover)
# Check distribution of data
turnover %>%
ggplot(aes(x = stag, color = factor(event),
fill = factor(event))) +
geom_histogram(aes(y = ..density..), alpha = 0.5) +
geom_density(alpha = 0.05) + scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") + theme_minimal() +
theme(legend.position = "top")
# Check for censored observations
n <- dim(turnover)[1]
cat((n - sum(turnover$event))/n * 100, "% of observations are censored")
# Separate variable based on type
NUM_COLS <- c("stag", "age", "extraversion", "independ", "selfcontrol", "anxiety", "novator")
CAT_COLS <- c("gender", "industry", "profession", "traffic", "coach", "head_gender", "greywage", "way")
# Transform CAT_COLS into categorical type
for (COL in CAT_COLS){
turnover[COL] <- turnover[COL] %>% unlist() %>% factor()}Number of missing values : 0Number of duplicates: 13
49.82079 % of observations are censoredResults
# Categorical covariates
# Selection of covariates (that are only discrete)
turnover.cat <- turnover %>%
select_if(is.factor) %>%
mutate(event = turnover$event)
turnover.cat %>% ggplot(aes(x = gender, color = factor(event),
fill = factor(event))) +
geom_bar(alpha = 0.5) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
turnover.cat %>% ggplot(aes(x = industry, color = factor(event),
fill = factor(event))) +
geom_bar(alpha = 0.5) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
turnover.cat %>% ggplot(aes(x = profession, color = factor(event),
fill = factor(event))) +
geom_bar(alpha = 0.5) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
turnover.cat %>% ggplot(aes(x = traffic, color = factor(event),
fill = factor(event))) +
geom_bar(alpha = 0.5) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
turnover.cat %>% ggplot(aes(x = coach, color = factor(event),
fill = factor(event))) +
geom_bar(alpha = 0.5) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
turnover.cat %>% ggplot(aes(x = head_gender, color = factor(event),
fill = factor(event))) +
geom_bar(alpha = 0.5) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
turnover.cat %>% ggplot(aes(x = greywage, color = factor(event),
fill = factor(event))) +
geom_bar(alpha = 0.5) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
turnover.cat %>% ggplot(aes(x = way, color = factor(event),
fill = factor(event))) +
geom_bar(alpha = 0.5) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
# Continuous covariates
turnover %>% ggplot(aes(x = stag, color = factor(event),
fill = factor(event))) +
geom_histogram(aes(y = ..density..), alpha = 0.5) +
geom_density(alpha = 0.05) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
turnover %>% ggplot(aes(x = event, color = factor(event),
fill = factor(event))) +
geom_histogram(aes(y = ..density..), alpha = 0.5) +
geom_density(alpha = 0.05) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
turnover %>% ggplot(aes(x = age, color = factor(event),
fill = factor(event))) +
geom_histogram(aes(y = ..density..), alpha = 0.5) +
geom_density(alpha = 0.05) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
turnover %>% ggplot(aes(x = extraversion, color = factor(event),
fill = factor(event))) +
geom_histogram(aes(y = ..density..), alpha = 0.5) +
geom_density(alpha = 0.05) + scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") + theme_minimal() +
theme(legend.position = "top")
turnover %>% ggplot(aes(x = independ, color = factor(event),
fill = factor(event))) +
geom_histogram(aes(y = ..density..), alpha = 0.5) +
geom_density(alpha = 0.05) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
turnover %>% ggplot(aes(x = selfcontrol, color = factor(event),
fill = factor(event))) +
geom_histogram(aes(y = ..density..), alpha = 0.5) +
geom_density(alpha = 0.05) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
turnover %>% ggplot(aes(x = anxiety, color = factor(event),
fill = factor(event))) +
geom_histogram(aes(y = ..density..), alpha = 0.5) +
geom_density(alpha = 0.05) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
turnover %>% ggplot(aes(x = novator, color = factor(event),
fill = factor(event))) +
geom_histogram(aes(y = ..density..), alpha = 0.5) +
geom_density(alpha = 0.05) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
# Other encoding of discrete variables
turnover.num <- turnover.cat[-9]
levels(turnover.num$gender) <- 1:length(levels(turnover.cat$gender))
levels(turnover.num$industry) <- 1:length(levels(turnover.cat$industry))
levels(turnover.num$profession) <- 1:length(levels(turnover.cat$profession))
levels(turnover.num$traffic) <- 1:length(levels(turnover.cat$traffic))
levels(turnover.num$coach) <- 1:length(levels(turnover.cat$coach))
levels(turnover.num$head_gender) <- 1:length(levels(turnover.cat$head_gender))
levels(turnover.num$greywage) <- 1:length(levels(turnover.cat$greywage))
levels(turnover.num$way) <- 1:length(levels(turnover.cat$way))
turnover.num <- as.data.frame(apply(turnover.num, 2, as.numeric))
turnover.num <- cbind(turnover.num, turnover %>%
select_if(is.numeric))
# Correlation matrix
corrplot(cor(turnover.num), col = brewer.pal(10, "BrBG"), method = "square", diag = FALSE)
Survival Plots
Now we will plot our categorical variables (gender, industry, profession, traffic, coach, head_gender, greywage, and way) using the Kaplan-Meier Estimation Curve. The Kaplan-Meier allows us to visualize how each category of survival probability differs from the others.
Exploratory Data Analysis
Kaplan-Meier Curve Stratified by Treatment Groups
# Estimate the median survival time among the three groups
survfit(Surv(time,status) ~ rx, data = colon_surv)Call: survfit(formula = Surv(time, status) ~ rx, data = colon_surv)
n events median 0.95LCL 0.95UCL
rx=Obs 315 168 2083 1656 2789
rx=Lev 310 161 2152 1540 NA
rx=Lev+5FU 304 123 NA 2725 NA# count the number of events after 2080 days, which is the median survival time among the observation group
tt <- colon_surv%>%filter(time > 2083)%>% group_by(rx)%>%summarise(ct = n(),
death = sum(status))
# Plot survival curve
fit <- survfit(Surv(time,status) ~ rx, data = colon_surv)
ggsurvplot(fit, data=colon_surv, risk.table = TRUE)
# Estimate the probability of surviving beyond 3000 days
summary(survfit(Surv(time, status) ~ rx, data = colon_surv), times = 3000)Call: survfit(formula = Surv(time, status) ~ rx, data = colon_surv)
rx=Obs
time n.risk n.event survival std.err lower 95% CI
3.00e+03 6.00e+00 1.68e+02 4.08e-01 3.97e-02 3.37e-01
upper 95% CI
4.94e-01
rx=Lev
time n.risk n.event survival std.err lower 95% CI
3.00e+03 4.00e+00 1.61e+02 3.92e-01 5.63e-02 2.96e-01
upper 95% CI
5.20e-01
rx=Lev+5FU
time n.risk n.event survival std.err lower 95% CI
3.00e+03 7.00e+00 1.23e+02 5.61e-01 3.43e-02 4.97e-01
upper 95% CI
6.32e-01 # compare significant difference in survival times between the three groups
survdiff(Surv(time, status)~ rx, data = colon_surv)Call:
survdiff(formula = Surv(time, status) ~ rx, data = colon_surv)
N Observed Expected (O-E)^2/E (O-E)^2/V
rx=Obs 315 168 148 2.58 3.85
rx=Lev 310 161 146 1.52 2.25
rx=Lev+5FU 304 123 157 7.55 11.62
Chisq= 11.7 on 2 degrees of freedom, p= 0.003 K-fold Cross Validation
#| echo: False
#| message: false
#| warning: false
#| include: False
library(survival)
library(caret)
# Fit the Cox model
cox_model <- coxph(Surv(time, status) ~ rx + age + surg + strata(obstruct) + strata(differentiation) + node4 + local_spread, data = df)
# Calculate the original c-index
c_index_original <- survConcordance(Surv(time, status) ~ predict(cox_model), data = df)
#cat("Original c-index:", c_index_original$concordance, "\n")
# Create a function for calculating c-index in each fold using survConcordance
cox_cindex <- function(train_data, test_data) {
fit <- coxph(Surv(time, status) ~ rx + age + surg + strata(obstruct) + strata(differentiation) + node4 + local_spread, data = train_data)
# Calculate concordance on test data
c_index <- survConcordance(Surv(test_data$time, test_data$status) ~ predict(fit, newdata = test_data), data = test_data)$concordance
return(c_index)
}
# Perform 5-fold cross-validation with stratification
K <- 5
folds <- createFolds(df$status, k = K, list = TRUE, returnTrain = TRUE)
cv_c_indices <- sapply(folds, function(train_indices) {
train_data <- df[train_indices, ]
test_data <- df[-train_indices, ]
cox_cindex(train_data, test_data) # use the concordance function inside cox_cindex
})
# Print cross-validated c-indicesOriginal c-index: 0.6544784 Mean cross-validated c-Index: 0.6420104 